Relativistic hydrodynamics in the presence of puncture black holes 



o 
o 

(N 
o 

o 

in 



o ; 



> 
m 

00 
O 

o 



X 



Joshua A. Faber,^'EI Thomas W. Baumgarte,^'[l| Zachariah B. Etienne/ Stuart L. Shapiro, ^'0 and Keisuke Taniguchi^ 

^Department of Physics, University of Illinois at Urbana- Champaign, Urbana, IL 61801 
^Department of Physics and Astronomy, Bowdoin College, Brunswick, ME O4OII 

(Dated: February 1, 2008) 

Many of the recent numerical simulations of binary black holes in vacuum adopt the moving 
puncture approach. This successful approach avoids the need to impose numerical excision of the 
black hole interior and is easy to implement. Here we wish to explore how well the same approach can 
be applied to moving black hole punctures in the presence of relativistic hydrodynamic matter. First, 
we evolve single black hole punctures in vacuum to calibrate our BSSN (Baumgarte-Shapiro-Shibata- 
Nakamura) implementation and to confirm that the numerical solution for the exterior spacetime is 
invariant to any "junk" (i.e., constraint-violating) initial data employed in the black hole interior. 
Then we focus on relativistic Bondi accretion onto a moving puncture Schwarzschild black hole as 
a numerical testbed for our high-resolution shock-capturing relativistic hydrodynamics scheme. We 
find that the hydrodynamical equations can be evolved successfully in the interior without imposing 
numerical excision. These results help motivate the adoption of the moving puncture approach to 
treat the binary black hole-neutron star problem using conformal thin-sandwich initial data. 



I. INTRODUCTION 

The evolution of spacetimes containing black holes 
(BHs) has been a long-standing focus of numerical rel- 
ativity. Among other astrophysical systems containing 
compact objects, BHs in binaries with compact object 
companions have long been considered promising sources 
of gravitational waves that can be detected by the new 
generation of ground-based gravitational wave interfer- 
ometers LIGO 1], TAMA _2], GEO Q and VIRGO [1], 
and by the proposed space-based detector LISA 0]. 

The recent association of short gamma-ray bursts with 
galaxies with extremely low star- format ion rates has vir- 
tually ruled out supernovae as the progenitors of short 
gamma-ray bursts, and instead favors either stellar-mass 
black hole-neutron star (BHNS) or neutron star-neutron 
star (NSNS) mergers (see p for a review). While mul- 
tiple aspects of the NSNS scenario have been studied in 
detail (see, e.g., 0, i, i, O, [O, El, El for simulations 
performed in general relativity), progress on BHNS merg- 
ers has been significantly slower. This is not surprising, 
since simulations of BHNS binaries combine the compu- 
tational difficulties associated with the BH singularity 
with those arising from shocks and other hydrodynamic 
phenomena associated with the neutron star matter. 

Recent advances in the numerical simulation of BHBH 
binaries (see El, El| as well as numerous follow-up 
papers) have overcome many of the difficulties associated 
with the BH singularity. In particular, the "moving punc- 
ture" approach adopted by [H, El does not require any 
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excision of the BH interior, is quite easy to implement 
into BSSN numerical relativity schemes [13, El, and has 
been used successfully in a large number of simulations. 
These simulations have treated the case of equal-mass, 
non-spinning binaries [l^, [13, HU, [H, [ll , and in some 
cases have also dealt with binaries with unequal mass 
0, HE HI, non-zero spins either aligned with the or- 
bital angular momentum [l^, [H, [H, [11 or at arbitrary 
orientations 31, 32, 33], or combinations of both [s^. [Sq 
(see also [36|, (STI |38|, |39| for alternative approaches to the 
BHBH problem). These puncture simulations have fur- 
nished results of great astrophysical interest, including 
the gravitational wave signals expected from such merg- 
ers and the kick velocities imposed on the merger rem- 
nant. 

In the moving puncture approach the presence of the 
singularity is largely ignored, which raises the question 
as to why it does not spoil the numerical evolution. This 
issue has been clarified by [iO, [4l[ , who analyzed the geo- 
metric structure of puncture solutions. Under the gauge 
conditions that are used in moving puncture simulations, 
the dynamical evolution of a single isolated Schwarzschild 
BH settles down to a time slicing that terminates at a 
limiting surface of finite areal radius. This slicing there- 
fore avoids the central curvature singularity, and covers 
only regular regions of the Schwarzschild geometry. In 
effect, the moving puncture approach provides a means 
of "excision-without-excision" (see also [H, [11 and Sec- 
tion HTC] below). An analytic expression describing the 
asymptotic, late-time solution has been found in [S| , and 
we use this solution to test the convergence behavior of 
our code in the presence of a puncture BH. 

The fact that time slices in the moving puncture ap- 
proach cover only regular regions of spacetime makes 
it an attractive approach for modeling BHNS binaries, 
since then the equations of relativistic hydrodynamics - 
or any other matter model - can be integrated together 
with the gravitational field equations without any need 
for excision. In fact, this is the approach adopted in the 
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only fully self-consistent dynamical simulations of orbit- 
ing BHNS binaries to date ^,\^ (see gUii, [H, [H for 
other preliminary relativistic simulations of BHNS bina- 
ries, and for discussion about and simulations of 
alternate configurations without the use of excision). 

In this paper we analyze in detail the hypothesis that 
relativistic hydrodynamics in the presence of puncture 
BHs does not require excision. As a test problem we 
study relativistic Bondi accretion onto a BH and compare 
with the analytic solution for stationary, spherical flow. 
We perform this test both for static and moving black 
holes (the latter at rest with respect to the asymptotic 
gas flow), and verify that the fluid behaves as expected. 

This paper forms a natural stepping-stone in our 
group's systematic effort to model BHNS binaries. Af- 
ter solving the initial value problem in general relativity 
[53I [5^ . [55I . [56j and performing preliminary dynamical 
simulations in conformal gravitation [47l . |48 l|, we have 
now assembled and tested the tools to model BHNS bi- 
naries self-consistently. In contrast to [i^ |4| , who con- 
structed BHNS initial data with a puncture method (see 
Section III Bl below), we plan to evolve our quasiequilib- 
riuni initial data, which we have constructed via the con- 
formal thin-sandwich method. 

Our forthcoming evolution of conformal thin-sandwich 
initial data raises another conceptional issue that we ad- 
dress in this paper. To construct quasiequilibrium BH 
initial data requires excising the interior and imposing 
suitable boundary conditions on the BH horizon [s^, ■ 
This approach provides data only in the black hole exte- 
rior, but the moving puncture approach for the dynam- 
ical evolution requires data in the BH interior as well. 
However, by definition of the BH horizon, data in the BH 
interior cannot affect the exterior spacetime. We demon- 
strate in Section IIII Al that we can indeed replace data 
in the BH interior with constraint-violating "junk" data, 
and still find an identical evolution in the exterior. In 
fact, the evolution settles down to the same asymptotic 
late-time solution in the BH interior as well, indepen- 
dently of the initial data. These findings thus suggest 
that moving puncture simulations of BHBH binaries can 
adopt conformal thin-sandwich initial data, which may 
be less inherently eccentric and thus be more accurate 
approximations to true quasi-equilibrium than the punc- 
ture initial data typically used in moving puncture binary 
simulations [59"|. This hypothesis is tested in detail in a 
forthcoming work [60]. 

This paper is organized as follows. In Sec.|TTl we sum- 
marize the puncture BH formalism, as well as our specific 
implementation of initial data, gauge conditions, matter 
evolution, and code diagnostics. In Sec. IIIIl we discuss 
results for a number of vacuum spacetime evolutions per- 
formed for both stationary and moving puncture BHs, 
including calculations for which we alter the initial data 
in the BH interior. In Sec. IIVI we present simulations 
of puncture BH spacetimes that contain hydrodynamic 
matter, and discuss how we implement a scheme that al- 
lows for stable and accurate hydrodynamical evolutions. 



We conclude in Sec. |V] with a brief discussion of our find- 
ings. 

II. NUMERICAL IMPLEMENTATION 

A. The 3+1 decomposition 

Throughout this paper we will cast the spacetime met- 
ric gab into the 3+1 ADM (Arnowitt-Deser-Misner) form 

ds^ = ~{a^ - P,f3')dt'^ + 2(3idtdx' + ^'^%dx'dx^, (1) 

where a is the lapse function, the shift vector, ip the 
conformal factor, and 7^ = ttio conformally re- 

lated spatial metric, defined so that its determinant 7 = 1 
in Cartesian coordinates. The extrinsic curvature Kij is 
defined by 

idt-C(3)-f^j=-2aK,j, (2) 

where £ denotes the Lie derivative. We note that two 
independent conformal rescalings of the tracefree part of 
the extrinsic curvature Aij are widely used in the litera- 
ture. In the context of initial data decompositions, this 
quantity is usually rescaled as 

A, ^ (^Ku - ^g^.K^ , (3) 

whereas BSSN-style evolution schemes typicahy rescale 
■Aij as 

4- ^ V^-' {k,j - l9^,K^ , (4) 

so that its indices may be raised and lowered with the 
conformally related spatial metric. Our code evolves the 
latter expression. 

We adopt the BSSN formulation [13, [11 for the dy- 
namical evolution of the gravitational fields (see also 
Eqs. (11) - (15) of [lH). In addition to the quantities 
7ij, (j) = In?/', Aij and K, this formulation utilizes the 
"conformal connection functions" = — 7*j as auxiliary 
quantities. Like many BSSN implementations, we en- 
force the vanishing trace of Aij and unit determinant of 
7jj at every timestep. 

B. Puncture initial data 

The basic idea of puncture initial data is to factor 
out from the spatial metric analytic terms that repre- 
sent the singular terms at a BH singularity, and treat 
only the remaining regular terms numerically |6^.[63l[6^. 
Specifically, under the assumption of conformal flat- 
ness and maximal slicing, the momentum constraint de- 
couples from the Hamiltonian constraint in the confor- 
mal transverse-traceless decomposition of Einstein's con- 
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straint equations (see, e.g., [6a, ISa] for reviews). More- 
over, the momentum constraint becomes linear and a so- 
lution representing a BH with momentum is given by 



(7'^-nV)F'=nfe], (5) 



where r and n* are the distance to and radial unit nor- 
mal vector from the BH, and the expression is scaled as 
in Eq. ([3]). This expression can then be inserted into the 
Hamiltonian constraint. Decomposing the conformal fac- 
tor as a sum of an analytic, singular part tps = 1 4- /2r 
and a term u that is regular everywhere. 



1 



M 

2r ' 



(6) 



where M. is a, constant, the Hamiltonian constraint be- 
comes a regular equation for u, 
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(7) 



For moving BHs, the solution for ip can be found analyt- 
ically up to second order in (see, e.g., (67| ) 

= 8M{M + 2rf [""('^^-Po^^^) + "2(0^2(a^)] (8a) 

M_ 



U2 



3787^"^ 



13167W''r 



3^3 



-^21567^^"* -t- 1536Alr^ + 240r^ 



-217W(X + 2r)^lni 



M 



M + 2r' 



(8c) 



where Po{fi) = 1 and P2(m) = 3/i^/2 — 1/2 are Legen- 
dre polynomials and /i = cos0. Clearly, for static BHs 
with = we find u = as expected, and recover a 
Schwarzschild t = const, slice in isotropic coordinates. 
We also note that the parameter A4 in the above equa- 
tions reduces to the BH ADM mass only for static BHs 
(see also Section Hi Gl below) . 



With suitable coordinate conditions (see Section [HD] be- 
low), this prescription leads to remarkably stable evolu- 
tions. This is somewhat surprising, since one might have 
expected the presence of singularities to spoil the numer- 
ical evolution. This issue has been clarified recently by 
pol , who analyzed the geometric structure of punc- 
ture solutions. 

For a single isolated BH at rest in the coordinate 
system, initial data representing a slice of constant 
Schwarzschild time in isotropic coordinates are given by 
Eq. ^ with u = 0. Isotropic coordinates do not pen- 
etrate the BH horizon, and instead cover two copies of 
the BH exterior that can be thought of as two different 
asymptotically flat universes connected by an Einstein- 
Rosen bridge. The conformal factor Eq. ([6]) diverges at 
r = 0; this point corresponds to the asymptotically flat 
end of the "other" universe, and hence is a relatively 
harmless coordinate singularity only. This initial slice 
does not encounter the BH spacetime singularity. 

With the gauge conditions typically used in moving 
puncture simulations, such initial data lead to an episode 
of "dynamical" evolution until the solutions settle down 
to a new equilibrium state. Specifically, in the adopted 
gauge the metric coefficients depend on time, and hence 
appear dynamical until they settle down. In the new 
equilibrium state the conformal factor is still singular at 
r = 0, but it now features a singularity instead of 

the previous 1/r coordinate singularity. Such behavior 
indicates that the point r = now represents a limiting 
surface of finite areal radius Parcai- The slice has thus 
disconnected from the other asymptotically flat end, and 
instead terminates on a surface of finite Paroai inside the 
BH's horizon. The numerical grid does not include the 
spacetime singularity at Parcai — 0, so that the point 
r = is again a coordinate singularity only. This helps 
explain the success of this numerical method, which, in 
short, may be thought of as "excision- without-excision" 
(see also [i^ for an explanation of this numerical behav- 
ior). 



C. Moving puncture evolutions 

The success of the puncture method for initial data 
suggests that a similar approach might be equally suc- 
cessful for dynamical evolutions. However, the original 
puncture method - namely, factoring out analytical sin- 
gular terms and evolving the remaining regular terms 
alone - did not achieve long-term stable evolutions in 
dynamical simulations (see, e.g., [H, [69|). The problem 
may be associated with the need for a coordinate system 
that leaves the puncture at a prescribed location in the 
numerical grid, given by the singularity in the analyti- 
cal function. The breakthrough in the recent dynamical 
puncture simulations is based on using a "movin g" p unc- 
ture in which no singular term is factored out [iSl . [l^ . 



This realization about puncture geometry suggests 
that it should be straightforward to extend the puncture 
approach to describe relativistic hydrodynamics in the 
presence of BHs. The numerical grid covers only regular 
regions of the spacetime, so the fluid cannot encounter 
any spacetime singularities. At r = 0, which corresponds 
to a sphere of finite areal radius Parcai inside the horizon, 
all field and fluid characteristics point inward to smaller 
areal radii. The interior of the sphere may therefore be 
disregarded - it cannot affect the exterior. As we will 
discuss in Section III El below, finite differencing around 
the point r — can lead to numerical error that compli- 
cates some of our expectations; however, this is a purely 
numerical issue that can be dealt with quite straightfor- 
wardly. 
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D. Gauge and boundary conditions 



E. Hydrodynamics 



The key to the success of a moving puncture evolution 
is the identification of suitable gauge conditions. Both 
[m and [iBl use an advective "l+log" shcing condition 
for the lapse 



dta - P'd.a = 2aK. 



(9) 



The advection term forces the singularity, located at the 
point where tjj oo and a — > 0, to move with coordinate 
velocity = dx^/dt = — We note that we evaluate 
all advective terms using "upwind" differencing, as 
described in [70l |. 

The gauge condition used in moving puncture evolu- 
tions is a hyperbolic "gamma-driver" condition for the 
shift, either 



dt(3' - (3^dj(3' = -B\ 



(10a) 



4- 

dtB'~l3^djB' = dtr ~ f3^ djf' - r,B' , (10b) 



dtP' = -^B\ (11a) 
dtB^ = dtr-r]B\ (lib) 

The former choice is referred to in [71| as the "shifting- 
shift" condition, and represents the choice Cs = Cb — 
Q — in the notation of [T^I, who studied the hy- 
perbolicity of the resulting evolution scheme. The lat- 
ter is referred to as the "non-shifting-shift" , and has 
Cs = Cb = C/ = 1- We note that [4g use a modified, 
first-order version of the non-shifting-shift with a differ- 
ent stabilization term. Throughout this work, we will use 
only the "non-shifting -shif t" , E qs. pTa|) and (fTTb)) . 

The term rj in Eqs. (jlObp and (jllbp . which has units of 
A4~^, is a damping term that has a non-trivial effect on 
the evolution of the puncture as it moves across the grid. 
In general, for larger values of 77, the BH has a smaller 
coordinate velocity but a larger coordinate radius, which 
allows for better numerical resolution of the BH horizon. 
However, larger values of rj were shown in [7l!| to produce 
significantly larger values of F*, indicating a stronger de- 
formation of the metric during the evolution. This is 
often disadvantageous, especially with regard to specify- 
ing boundary conditions. We use 77 — 0.25/Al for the 
static BHs in Sec. IHI Al rj = 0.2/M. for matter evolu- 
tions in Sec. IIVI but rj ~ 2.0/A4 for moving black holes 
in Secs. lIHB ll and llHB 2l so that we can directly compare 
with the results of [71]. 

For our runs with relatively close outer boundaries de- 
scribed in Sec. IHI Al we use Robin-type boundary condi- 
tions for our metric and field variables, assuming a l/r 
falloff behavior away from their asymptotic values. For 
larger runs we employ outgoing wavelike boundary con- 
ditions, based on the local speed of light and again as- 
suming a l/r falloff. The conditions are slightly modified 
when we use a "fisheye" coordinate scheme, as we discuss 
in Appendix [XI 



Our hydrodynamics scheme is essentially identical to 
the grid-based, fully general relativistic, high resolution 
shock capturing scheme described in [73|, although here 
we do not include any electromagnetic terms in our evo- 
lution calculations. This scheme for both matter and 
field evolution is based on second-order finite differenc- 
ing, which we have now embedded within the Cactus grid 
hierarchy f?^. 

We assume a stress energy tensor in the form T'^'^ — 
Pohu'^u'^+Pg'^^ where po is the rest density, h the specific 
enthalpy, and u^^ the 4-velocity of the fiuid. We will also 
assume a gamma-law equation of state 



P= (F-l)poe, 



(12) 



where e is the specific internal energy. 

In our numerical code, we evolve the "conserved" quan- 
tities yO*, Si, and f, defined by Eqs. (36) and (37) of 
[73| . These conserved quantities are defined in terms of 
the "primitive" variables p, w*, and P, according to their 
Eqs. (41) ~ (43). The transformation from the conserved 
quantities back to the primitive set requires an iteration 
(see Eqs. (59) - (62) of [zl). 

We find that this iteration is very accurate everywhere 
except in the region immediately surrounding the punc- 
ture, where accumulated errors in the finite differencing 
can lead to unphysical values of the primitive variables. 
To overcome this difficulty, we note that the quantity f, 
which serves as the energy variable in the conserved set 
of variables, should always remain positive. In fact, it 
has to remain greater than w — p*, where w = avP p^:, in 
order to satisfy h>l (see Eq. (58) in [TSj]) , but positiv- 
ity is sufficient for our purposes. To enforce this, we set 
f to a minimum value, typically f — whenever it 

would otherwise be negative. Furthermore, from (58) in 
[7^ we see that in the limit /i ^ 1 we have -u; ^ f -t- , 
or equivalently 



\S\^ = ^''S,S,^t{t + 2p,). 



(13) 



It can be shown that this limit on \S\'^ is an upper limit. 
When numerical error leads to a value of \S\'^ that is 
larger than the allowed value, we enforce this condition 
by reducing the magnitude of Si so that 15^ is reduced 
to 0.98 X f (f -I- 2/9*), leaving the ratios of the individual 
vector components unchanged. 

We again point out that the above "fixes" are employed 
only in the immediate vicinity of the puncture, inside the 
BH horizon. The need for these fixes arises from poor res- 
olution around the puncture and the resulting large finite 
differencing error. We find that increasing the grid reso- 
lution decreases the size of the region in which the fixes 
are needed. Other than dealing with this break-down in 
the iteration between conserved and primitive fluid vari- 
ables, the equations of relativistic hydrodynamics can be 
integrated without need for excision in the presence of 
puncture black holes. 
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F. Coordinates 

We perform simulations both in 2+1 (i.e., axisymme- 
try) and in 3+1, using the HRSC, BSSN code of [zl, 
which accommodates both cases. Our 2+1 simulations 
adopt the "cartoon" method for the fields, as described 
in [zl. 

Finite difference simulations with fixed spatial resolu- 
tion in three spatial dimensions are very intensive com- 
putationally, since they simultaneously require a fine grid 
resolution close to the BH, and boundary conditions im- 
posed at a sufficiently large separation. We have imple- 
mented so-called "fisheye" coordinates to solve the prob- 
lem of dynamic range (see [l5| ) . In Appendix [X] we re- 
view the coordinate transformation and its effect on the 
various terms that appear in our numerical scheme, es- 
pecially our shift evolution equation. 

G. Diagnostics 



advective term j3^dia. Most commonly, moving punc- 
ture solutions employ condition Eq. ^ with the advec- 
tive term; we will refer to this condition as "advective 
1+log" slicing. The geometry of the resulting spacetime 
has been analyzed in In [4l|, the authors consid- 

ered dropping the advective term, which results in the 
"non-advective 1+log" slicing 

dta = -2aK. (18) 

As it does for the advective 1+log slicing, the dynami- 
cal evolution using the modified slicing condition quickly 
settles down to a new equilibrium, but in this case this 
time-independent solution must evidently be maximally 
sliced with K = 0. As it turns out, this maximally sliced 
asymptotic solution can be found analytically (see 44]). 
In Section 11101 we therefore perform tests comparing 
with the analytic solution, before considering the more 
common advective 1+log slicing in Section UlIBI 



We monitor several global quantities computed from 
surface integrals at large separations. These include the 
ADM mass 

Madm = ^ / (rV8 - Y'd,ij)dj:,, (i4) 

and the linear momentum 

Pj^^j>{K]~5]K)d^,, (15) 

where dSi = {x'^ /r)ip^r'^ sinOdOdcf) for a spherical surface 
at fixed radius. 

We also compute the irreducible mass of the BH from 

Mi„ = y^A/Wn. (16) 

Instead of computing the proper area A of the BH's event 
horizon, we approximate this area as that of the apparent 
horizon. 

For the approximate boosted BH solutions of Section 
III B I the ADM mass and irreducible mass are related to 
the parameter M, to leading order in P% by 

5 p2 I p2 

Madm = M + - —, M-„,= M + - (17) 

(see [zl]). Throughout this paper we express dimensional 
quantities in units of M. 

III. VACUUM TESTS 

To test the vacuum sector of our puncture code, we 
have performed a suite of tests for stationary and mov- 
ing puncture BHs, both in 2+1 and 3+1. In these tests 
we adopt two slightly different slicing conditions, namely 
the 1+log condition Eq. ([5]), both with and without the 



A. Non-advective 1+log slicing 

As demonstrated in [4l| . dynamical puncture evo- 
lutions with the non-advective 1+log slicing condition 
Eq. (fT8|) settle down to a maximally sliced, time- 
independent solution, given analytically by Eqs. (11)-(15) 
of [AA\ . We test our code by comparing the late-time solu- 
tion of dynamical puncture evolutions with the analytical 
solution. In these tests we consider three different types 
of initial data, which we summarize in Table [H In the 
first set of runs, labeled by "BN" in TableHl we adopt the 
analytical solution of [4J] itself as initial data. The sec- 
ond set of runs, labeled by "ISO", adopts as initial data 
& t — const, time slice of the Schwarzschild solution in 
isotropic coordinates, given by Eq. ([6]) for u = 0, together 
with a = ip^'^. The third set of initial data is identical to 
the second set in the BH exterior, but now we replace the 
interior data with some constraint- violating "junk" . We 
adopt these data to test whether we can use initial data 
that only provide valid exterior data in puncture evolu- 
tions, which also require initial data to fill the arrays in 
the interior. We consider three different choices for the 
interior solution. In the first choice, labeled "Horizon 
Junk" or "HJ" , we set the values of the conformal factor 
and lapse function throughout the interior of the BH to 
their values on the horizon, i.e. tp = 2 and a = 1/4 for 
r < Th = Ai /2. For the second choice, labeled "Interior 
Junk" or "IJ" , we again set ip and a to constant values, 
but this time only inside rh/2^ so ip = 3 and a = 1/9 for 
r < r/i/2, using the isotropic solution everywhere outside 
this. Finally, in our third choice, labeled "Smooth Junk" 
or "SJ" , we construct an even, fourth-order polynomial 
in the BH interior that joins on to the exterior conformal 
factor in such a way that the function and its first two 
derivatives are continuous on the horizon. Specifically, 
we choose -0 = 2.875 - 5{r/M)^ + 6{r/M)-^ in the BH 
interior, and a = . 
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TABLE I; Summary of our stationary vacuum puncture simu- 
lations discussed in Sec. lIII A1 "Baumgarte & Naculich" refers 
to the analytical solution of ,4^]; our different "junk" initial 
data are described in the text. 



Name 


Initial Data 


Resolution 


Grid 


2+1; no fisheye; ji* 


< IbM 


BN2a 


Baumgarte & Naculich 


X/32 


480 X 480 


BN2b 


Baumgarte & Naculich 


M/24 


360 X 360 


BN2c 


Baumgarte & Naculich 


M/IG 


240 X 240 


BN2d 


Baumgarte & Naculich 


M/S 


120 X 120 


HJ2a 


Isotropic "Horizon Junk" 


X/32 


480 X 480 


IS02a 


Isotropic Schwarzschild 


X/32 


480 X 480 




3+1; fisheye, \x'\ < QM - 


\x'\ < 15M 


IS03 


Isotropic Schwarzschild 


M/16 


192^ X 96 


IJ3 


Isotropic "Interior Junk" 


M/16 


192^ X 96 


HJ3 


Isotropic "Horizon Junk" 


M/16 


192^ X 96 


SJ3 


Isotropic "Smooth Junk" 


M/16 


192^ X 96 



All simulations described in this Section were per- 
formed using equatorial symmetry. For 3+1 runs we 
used fisheye coordinates with parameters oq = 1, ai = 4, 
i?i = 3, si = 1.0. While we always employ a cubic coor- 
dinate grid, a cubic grid in fisheye coordinates does not 
correspond to a cubic grid in physical coordinates. Thus, 
the boundaries we quote in physical coordinates only ap- 
ply at the coordinate axes, and lie further out at other 
angles. When we quote a numerical resolution for a given 
run, it refers to the central region of the fisheye transition, 
for which the physical and fisheye coordinates essentially 
overlap. The numerical resolution with respect to phys- 
ical coordinates will be more coarse at larger radii, by 
construction. In all calculations described in this sec- 
tion, we set 7] — 0.25/ A4 in the shift evolution equation. 
All evolutions described in this Section were terminated 
at tp = 3007W. 

We first compare the "BN2" simulations at four dif- 
ferent grid resolutions to establish second-order conver- 
gence of our code. In the top panel of Fig. [U we show 
the absolute deviations between our numerical runs at 
t = IQM and the exact solution for run BN2a at reso- 
lution M/32 (solid), BN2b at 7W/24 (dotted), BN2c at 
M/16 (dashed), and BN2d at M/8 (dot-dashed), focus- 
ing on the inner region of the grid, where deviations are 
most apparent, and which, ai t — 10A4, is causally dis- 
connected from the outer boundary at 15 A4. To confirm 
the expected second-order convergence we show the same 
quantities in the bottom panel rescaled by factors of 16, 
9, 4, and 1, respectively. 

We now compare the late-time solution as it emerges 
from our three different types of initial data. In Fig. [5] we 
show results for runs BN2a, IS02a and IIJ2a, which all 
have the same grid resolution. Evidently, all late-time so- 
lutions agree very well with the analytical solutions. The 
deviation from the analytical solution is at most about 
1% for all three runs, and is caused by finite difference 
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FIG. 1: Deviations from the exact solution for the conformal 
factor Ai/) for the runs BN2 (see Table HI. The top panel 
shows the deviations themselves, while in the lower panel the 
deviations are appropriately rescaled to demonstrate second- 
order convergence. 



error as well as the outer boundary. To within these 
numerical errors, all three late-time solutions agree with 
each other, even if we replace the interior initial data 
with junk that does not even satisfy the constraints. It 
is particularly noticeable that even in the BH interior no 
"memory" of this junk remains at this late time, and that 
the solution approaches the analytically known late-time 
solution throughout. 

While it is not surprising that data in the BH interior 
does not affect the exterior solution, it is reassuring to see 
that this holds true even in a numerical simulation. To 
further illustrate this point we performed the 3+1 simula- 
tions IS03, IIJ3, IJ3, and SJ3, and plot the ADM mass as 
a function of time in the top panel Fig. [31 We find excel- 
lent agreement with the stationary solution at late times 
for these runs, indicating that information contained in 
the BH interior does not affect the exterior solution as 
we evolve. In quantitative terms, applying junk either 
within the horizon or smoothly at the horizon results in 
a change in the ADM mass of 0.05%, as we show in the 
bottom panel of the Figure. We note that reflections off 
the outer boundary are at least partially responsible for 
the discrepancy, given the time at which they appear. To 
confirm the proper behavior of our code when we have 
constraint-violating junk in the BH interior, we show in 
Fig. m the convergence of the ADM mass for stationary 
puncture runs using smooth junk, performed at numer- 
ical resolutions of A^/20, A^/16, and M/12. The lower 
panel suggests self-consistent second-order convergence 
at all times. At both early and late times, we have also 
verified that the code converges at second-order to the 
analytic solution (Madm = 1); at intermediate times we 
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300 



FIG. 2: The top panel shows profiles of the conformal factor 
V) at f = 300X for for runs BN2a, IS02a and HJ2a (see 
Table All simulations are performed with a grid spacing of 
A^/32, using the gauge condition Eq. (I18|l . for which the late- 
time solution is known exactly. In the bottom panel we show 
the absolute deviations between our numerical results and the 
exact solution, finding that they are small throughout. The 
arrows mark the location of the apparent horizon. 



FIG. 3: The ADM mass (Hi} for the 3+1 runs using fish- 
eye and M/IQ numerical resolution in the central region (top 
panel). We show results for initial data IS03, IJ3, HJ3 and 
SJ3 (see Table|T]|. In the bottom panel, we plot the differences 
between the ADM mass from the junk runs and that of run 
IS03. We see a deviation of .05% between run IS03 and runs 
IJ3 and SJ3 over the course of the run, indicating that the 
evolution is insensitive to the initial data in the BH interior. 



find a brief departure from this behavior that is likely 
caused by reflections off the outer boundary. Even for 
the simulation HJ3, where the "junk" joins the exterior 
with discontinuous derivatives on the horizon, the ADM 
mass deviates from the IS03 simulations by only 1% over 
the course of the evolution. These simulations indicate 
that the evolution is relatively insensitive to the details 
of the initial data in the BH interior, even when the dif- 
ferencing stencil in the exterior of the BH overlaps non- 
differentiable regions. Still, our results clearly suggest 
that only smoothly extrapolated data should be used in 
dynamical evolution calculations. 

Our "junk" tests suggest a very simple solution to the 
following conceptional issue. Many sets of initial data de- 
scribing compact binaries are constructed with the con- 
formal thin-sandwich decomposition. For the construc- 
tion of BHs, these equations are supplemented with equi- 
librium boundary conditions that are imposed on the 
BHs' horizons. As a consequence, these data describe 
only the exterior geometry, and do not provide any ini- 
tial data for the BH interior. Evolution calculations that 
employ the moving puncture approach do not excise the 
BHs, and hence require initial data that extend into the 
BH interior. Our tests suggest that we can neverthe- 
less use excised initial data, e.g., conformal thin-sandwich 
data, and simply fill the BH interior with some arbitrary 
but sufficiently smooth "junk" . This is the approach that 
we plan to adopt for simulations of mixed BHNS binaries, 
enabling us to use our quasi-equilibrium models [ssl . [56j 
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FIG. 4: The top panel shows the ADM mass measured at a 
radius r = \2.\M, for runs with smooth junk as initial data 
and numerical resolutions of M /20, Al/16, and A^/12. In the 
bottom panel, the differences between the runs are rescaled 
to demonstrate second-order convergence. 



as initial data. 
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TABLE II: Vacuum evolution runs discussed in Sees. IIII B II 
and HUB 21 



Name Initial Data Grid 



Resolution tp 



2+1 w/Fisheye; \x^ish\ < 15M, \x'\ < 136X 



SP2a 


Stat. Punc. 


360 X 360 


7W/24 


IQQM 


SP2b 


Stat. Punc. 


240 X 240 


M/16 


lOQM 


SP2c 


Stat. Punc. 


120 X 120 


M/8 


IQOM 


3+1 w/Fisheye; < 12M, \x' \ < 88M 


SPSa 


Stat. Punc. 


240^ X 120 


M/IO 


100^1 


SP3b 


Stat. Punc. 


192^ X 96 


M/8 


lOOX 


SPSc 


Stat. Punc. 


144^ X 72 


M/6 


IQQM 


2+1 w/o Fisheye; = 0.5, \x'\ < 70M 


MP2a" 


Mov. Punc. 


1760 X 3520 


M/32 


bOM 


MP2b 


Mov. Punc. 


1680 X 3360 


M/24 


50X 


MP2c 


Mov. Punc. 


1120 X 2240 


M/16 


507W 


MP2d 


Mov. Punc. 


560 X 1120 


M/8 


507W 


3+1 w/Fisheye; 


= 0.5, 


< UM, \x' 


1 < 12QM 


MP3a 


Mov. Punc. 


336^ X 168 


M/V2 


^QM 


MP3b 


Mov. Punc. 


280^ X 140 


X/10 


507W 


MP3c 


Mov. Punc. 


224^ X 112 


M/8 


50A^ 


MP3d 


Mov. Punc. 


168^ X 84 


M/6 


50X 



Tor this run, < 55M 



B. Advective 1+log slicing 



The lapse evolution equation used in Sec. IIII Al above 
is useful since it leads to an analytically known solu- 
tion for a single BH. However, most moving puncture 
simulations of binaries employ the advective 1+log slic- 
ing Eq. ^ which allows the punctures to move through 
the numerical grid. We now study evolutions with this 
advective 1+log slicing, both for stationary and moving 
BHs. Our simulations are summarized in Table HIl where 
"SP" refers to stationary puncture solutions, described 
m more detail in Section HUB H and "MP" stands for 
moving puncture solutions, described in Sec. IIII B 21 For 
purposes of comparison comparison with the results of 
[nl (see Sec. HUB 21) we set r? = 2.Q/M in Eq. (fTTbll for 
all of the results discussed in this Section. 



our 2+1 and 3+1 simulations to study the convergence 
behavior of our code as indicated in Table |TT1 

In the top panel of Fig. [51 we show the ADM mass for 
our highest resolution 3+1 calculation, run SP3a, mea- 
sured at fint = 5.8, 6.8, 7.9, and 9.0. Similarly, in the 
bottom panel we show the ADM mass for our highest 
resolution axisymmetric calculation, run SP2a. We see 
two phenomena: first, at a time roughly corresponding 
to t = Tint , we see a small, temporary glitch in the ADM 
mass for each value of rjnt, followed by a slow deviation 
from the exact value. This indicates the passage of junk 
gravitational radiation, present in the initial data only 
because of numerical errors associated with discretiza- 
tion, through the surface, and is generic to calculations 
like these. Second, we note that our results converge as 
we move the integration surface outward, as we would 
expect. This is true both before and after the passage 
of the junk radiation through the surface. For the out- 
ermost surface, we see a variation in the ADM mass of 
no more that half a percent over an integration time of 
lOOAl. In general, we find that the late time deviations 
from the exact ADM mass scale like r^J^ through most 
of our grid. 
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1. Stationary black holes 

For our stationary puncture simulations we choose a 
double fisheye coordinate system with parameters ao = 1, 
ai = 4, 02 = 16, Ri = A.bM, R2 = 7.5A^, and si = S2 = 
1.5A^ (see Appendix]^. Our initial data are t = const, 
slices of the Schwarzschild metric expressed in isotropic 
coordinates, like the "ISO" data of Table HI We adopt 
the gauge conditions Eqs. (jllap and (jAlSp and evolve to 
a time tp = 100 Ai. By this time the solution has settled 
down into the late-time equilibrium solution, and further 
evolution would lead to only very small changes in the 
fields. We use three different grid resolutions for both 



FIG. 5: ADM mass versus time for stationary puncture evo- 
lutions, measured using Eq. (|14p at surfaces of different radii. 
In the top panel, we show results for surfaces at constant fish- 
eye radius f = 5.8 (dot-dashed), 6.8 (dashed), 7.9 (dotted), 
and 9.0 (solid) in 3+1 run SP3a, which correspond to phys- 
ical radii rmt = IIM, 17M, 27 M, and 42A^, respectively. 
In the bottom panel, the integration surfaces are placed at 
f = 5.8, 7.2, 8.5, and 9.9, corresponding to physical radii 
Tint = 11. IM, 19.9M, 35. 2A^, and 55. 2A^ in 2+1 axisym- 
metric run SP2a. For the outermost surfaces, we see devia- 
tions of less than half a percent over the course of the entire 
evolution in each case. 

To confirm convergence, we calculate the difference of 
the ADM mass from its analytical value as a function 
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of time for all the stationary puncture runs described in 
Tabic im In the 3+1 case, we take data from the integra- 
tion surface placed at nnt = 41.6, whereas for the 2+1 
case we choose the surface with nnt = 55.2. In Fig. [51 
we see the proper behavior for both the 3+1 and 2+1 
runs (top and bottom panels, respectively), confirming 
that our field evolution is indeed convergent to second 
order in the grid spacing. In the upper sub-panels, we 
show the expected convergence against the exact solu- 
tion. In both the 3+1 and 2+1 simulations, numeri- 
cal errors inherent in the finite-differenced initial data 
reach the integration surfaces at a time t « r-mf During 
the passage of oscillation, we still see second-order con- 
vergence in the differences between runs, plotted in the 
bottom sub-panels of each figure. Note that in the top 
panel, we rescale the higher resolution pair by a factor 
(6^2 _ g-2)/(-g-2 _ ^ 2.16, whereas in the bottom 

panel the rescaling factor is (8^^ — 16^^)/(16~^ — 24~^) = 
5.4. 
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FIG. 6: Differences in exact and numerical values for the 
ADM mass for our SP simulations of different numerical res- 
olution as a function of time (see Table |ll)| . In the top panel 
we show results for 3+1 simulations, with the ADM mass 
computed at rint = 41. 6A^. In the upper subpanel, results 
for different resolutions are properly rescaled to demonstrate 
second-order convergence against the exact solution. In the 
lower subpanel, we compare the differences between pairs of 
runs to demonstrate second-order convergence even during 
the passage of a spurious oscillation caused by initial errors 
through the data at rint • The bottom panel shows similar 
results for 2+1 simulations, with the ADM mass computed at 
Vint = 55.2M. Conventions are as above. 



Finally, in Fig. [71 we show the irreducible mass (see 
Eq. ((16])) of the BH as a function of time for our 3+1 
(top panel) and 2+1 (bottom panel) runs, for all three 
numerical resolutions in both cases. To determine the 
location and parameters describing the apparent horizon 



we use the Cactus thorn ahf inder [77| . We see that our 
results improve with increased resolution for both the 
2+1 and 3+1 simulations, as shown in the upper sub- 
panels. For our highest resolutions in the two cases, the 
deviation from the exact value are smaller than 0.5% and 
2%, respectively. In the lower subpanels, we show differ- 
ences between results for varying numerical resolutions, 
demonstrating second-order convergence until the hori- 
zon finder begins to show errors at late times for our 
lowest-resolution runs. 
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FIG. 7: Irreducible mass as a function of time for our SP runs 
(upper subpanels; see Table HI)) . For our highest resolution 
runs in 3+1 and 2+1, we see deviations of at most about 
2% and 0.5%, respectively. In the lower subpanels, we show 
the rescaled differences between pairs of runs, demonstrating 
second-order convergence. 



2. Moving Black Holes 

We now turn to moving BH solutions. In particular, we 
consider the parameters discussed in detail in ■71!], a BH 
with linear momentum = traveling along the 

X-axis, starting from an initial location xq — —3M. For 
our initial data, we calculate the extrinsic curvature and 
conformal factor from Eqs. (O, ([6]), and fSaj) - ([He]). We 
set the lapse initially to a{t = 0) = ip~'^{t = 0), and the 
shift to zero. The spatial metric is initially conformally 
flat. According to ^T7\ . the ADM mass of this config- 
uration is approximately AfADM — 1.15625A1, with the 
leading-order error term appearing at order P^. 

We use unigrid simulations for easier comparison with 
the results of [7l| , since the transformation to fishcyc co- 
ordinates would introduce new terms into the shift equa- 
tions (namely Eq. (|A15p instead of Eq. (|llbp ). Unigrid 
simulations are impractical in 3+1, so we perform ax- 
isymmetric simulations instead. We choose a numerical 
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grid that extends twice as far along the symmetry axis as 
it does radially. We note that while runs MP2b, MP2c, 
and MP2d, with numerical spacings Al/24, A4/16 and 
M /8, respectively, use a grid that extends to r = ±70A1 
along the symmetry axis (along which the BH moves) 
and r = 70A^ along the radial axis, our highest reso- 
lution calculation, run MP2a with spacing A4/32, only 
extends to r = 55A^ in both directions. We also perform 
3+1 simulations of moving punctures, which do use fish- 
eye coordinates. The parameters for these coordinates 
are the same as in Sec. IIII B 1[ except that we choose 
i?i = 6A^ and i?2 = 9A1 here to create a larger cen- 
tral region through which the BH travels. We evolve our 
moving puncture configurations until tp = 50A1, dur- 
ing which time the BH remains inside the central fisheye 
region. 

We show the ADM mass measured at various radii for 
our highest resolution 3-f 1 and 2-t-l simulations, MP3a 
and MP2a, in the top and bottom panels of Fig. [8] We 
again find better agreement between the measured ADM 
mass and the exact value when we evaluate Eq. at 
larger radii. 

Conservation of the ADM mass for moving punctures 
is comparable to that for stationary punctures, at ap- 
proximately 1% over a duration t = 50A^. Some of the 
initial deviations for the ADM mass in the 3-1-1 case, vis- 
ible especially for the innermost surface, are a result of 
small numerical errors associated with converting back 
into the physical value F* through a coordinate transfor- 
mation (see Appendix \K\ . 

As before, we confirm that the ADM mass converges 
to second order. In the top panel of Fig. [51 we show the 
ADM mass measured at a surface of radius rint = 62A1 
as a function of time for our four 3+1 runs performed 
using different numerical resolutions. As can be seen 
in Fig. [21 the values for the ADM mass slowly drift 
to larger values. We speculate that this is a result 
of the constraints equations being solved only to order 
P^, which introduces a small but non-zero error. In- 
stead of testing convergence to the approximate value 
of the ADM mass, we perform a self-consistence con- 
vergence test by by computing the difference in the 
ADM mass between pairs of runs in the bottom panel 
of the figure, scaling the results appropriately in each 
case. In terms of the difference between lowest resolu- 
tion pair, we scale the highest resolution pair by a factor 
(6-2 - 8-2)/(10-2 - 12-2) ^ 3 97 ^Yie medium reso- 
lution pair by a factor (6"^ _ ^^2y^^^-2 _ ^^-2^ ^ 2.I6. 
Our findings suggest second-order convergence through- 
out the evolution, even though at times after t=30M the 
analysis is complicated by small-amplitude numerical er- 
rors arising from discretization across fisheye transition 
regions. 

Combining the results from Figs. [51 and [HI we can ex- 
trapolate our results to both asymptotic radii and infi- 
nite numerical precision. With respect to the radius, we 
take our results at t = measured at r = 43A4 and 
62A^ for run MP3a, and assume leading-order 1/r falloff 



S 1.3 











R, ,-&ZM. 




r 


Ri„-4ait 






- R ,-2ait 






- R , -16.it 






ran MP3a (ii/12l - ' 1 _ ^ _ , - 




















R ;-5ait 






R i-44ii 






- R, ,-27it 






- R j-aoit 






ran MP2a (M./32) 

















10 20 30 40 50 



t/M, 

FIG. 8: ADM mass versus time for moving puncture 3+1 
evolution MP3a (see Table [IT]) , measured using Eq. ([14} at 
surfaces of different radii. In the top panel, we show results 
for surfaces at constant fisheye radius f = 8.0 (dot-dashed), 
9.2 (dashed), 10.5 (dotted), and 11.8 (solid), which corre- 
spond to physical radii nnt = 16 A^, 43 A^, and 62 A^, 
respectively. In the bottom panel, the integration surfaces are 
placed at nnt = 20X, 27 M, UM, and 53X for 2+1 simu- 
lation MP2a. For the outermost surfaces, we see deviations 
of less than a percent over the course of the entire evolution 
in each case, just as we did for the stationary puncture runs 
shown in Fig. [5] 



behavior in the measured ADM mass. With regard to 
numerical resolution, we Richardson extrapolate using 
runs MP3a and MP3b measured at r = 62A4. Com- 
bined, we find an extrapolated value for the ADM mass 
of A^ADM = 1.15615, which falls within 10-^ of the an- 
alytic value computed using Eq. ([14]) for these approxi- 
mate initial data. 

For our 3+1 simulations we also compute the linear 
momentum from Eq. (jl5p . In the top panel of Fig. 1101 
we show the linear momentum as a function of time, cal- 
culated at surfaces placed at radii nnt = 26.8M. In 
general, we see excellent agreement with the exact ana- 
lytic value up until t ~ rint, for the reasons noted above. 
Recall that we adopt initial data for moving BHs that 
are analytic but approximate, since they solve the con- 
straint equations only to order P^. The resulting error 
represents "junk" that propagates from the strong-field 
region outwards, and reaching the integration surface r^nt 
at approximately t ~ r^nt . Eq. (jlSp for the linear momen- 
tum also assumes that the coordinate vector pointing in 
the direction of the momentum is a Killing vector of the 
asymptotic, conformally related metric (see Appendix A 
of [2a]); this assumption breaks down once the "junk" 
reaches rint . Accordingly, evaluating Eq. (|15p leads to an 
error in the linear momentum, as is evident from Fig. 1101 
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FIG. 9: ADM mass versus time for moving puncture evolu- 
tions MPS (see Table HT]), for runs with different numerical 
resolutions (top panel). In each case, the integration sur- 
face was placed at r = 62M. In the bottom panel, we show 
the differences in the measured ADM mass between pairs of 
"neighboring" resolution, scaled to reflect the second-order 
convergence. 



FIG. 10: The top panel shows the linear momentum Px 
(Eq. (fT5)l ) evaluated at Tint = 26A4 for our simulations MPS. 
In the bottom panel, we show the difference between the val- 
ues computed with different numerical resolutions, scaling the 
differences between the more accurate runs by factors of S.97 
and 2.16, as in Fig. [§1 to highlight the second-order conver- 
gence of our code. 



In the bottom panel of the figure, we show the conver- 
gence behavior of the linear momentum. As for the ADM 
mass we find larger error terms once numerical errors 
reach the integration surface rjnt • 

As a further test we compare our numerical results for 
the lapse function a in runs MP2a and MP2c with those 
obtained by the numerical relativity group at NASA's 
Goddard Space Flight Center ([zH, hereafter GSFC). In 
Fig. [m we graph a along the trajectory of the puncture at 
t = AOAi ■ The simulations of GSFC was performed using 
a numerical resolution of A4/16, but features a fourth- 
order accurate differencing scheme. We nevertheless find 
excellent agreement between our results. 

Extending our tests to the 3-1-1 moving puncture cases, 
we show in the top panel of Fig. [T^] the lapse function 
along the axis on which the BH travels, at i = 20A4. 
In the bottom panel, we show the convergence of the 
lapse function by taking the difference between the runs 
and rescaling to second-order, finding good agreement 
throughout. The slight discrepancy near the puncture 
position is due to the fact that there is a sharp trough in 
the lapse function there, but the position of the puncture 
itself falls at slightly different coordinate positions for the 
three runs. Indeed, since the shift vector, like all field 
quantities, is second-order convergent in the numerical 
resolution, so too is the speed at which the puncture 
moves, since dx^/dt = — /3'(x*). 

Turning our attention to the position of the puncture, 
we show in Fig.[T3]the position of the BH puncture versus 
time in our three axisymmetric calculations performed at 
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FIG. 11: Lapse function a shown at T = 40 A4 along the path 
traveled by the moving puncture, for our axisymmetric 2-1-1 
runs with numerical resolutions of A4 /32 (run MP2a; dashed) 
and A4/16 (run MP2c; dotted), along with the numerical re- 
sults of GSFC (solid) as discussed in [7l[. The code used for 
the latter has numerical resolution M/16, but uses fourth- 
order differencing, whereas we use a second order scheme. 



different numerical resolutions, along with the Richard- 
son extrapolation value. We see good agreement, and 
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FIG. 12: The top panel shows the lapse function a measured 
on the X-axis, along which the puncture travels, at t = 20A^. 
In the bottom panel, we show the difference between the val- 
ues computed with different numerical resolutions, rescaling 
the differences as in Fig. [O] to highlight the second-order con- 
vergence of our code. 

note that numerical errors associated with coarser reso- 
lutions slow down the BH away from the proper asymp- 
totic velocity. In the bottom panel of the figure, we show 
the difference in position versus time for our 2-f 1 runs, 
again suggesting second-order convergence. 

IV. MATTER TESTS 

In this Section we describe relativistic hydrodynamics 
simulations in the presence of puncture BHs. We are par- 
ticularly interested in testing how the accretion of matter 
onto BHs can be simulated within the moving puncture 
approach. The exact Bondi solution for accretion onto a 
static Schwarzschild BH provides a perfect test bed for 
these purposes. A detailed discussion of the relativistic 
Bondi solution can be found in Appendix G of ^79]; we 
summarize all relevant expressions in Appendix [B) In 
Section flV Al we test our code's capability of simulating 
the accretion onto a static BH; in Section IIVBI we treat 
the identical problem, but viewed in a frame in which 
the BH is moving and represented by a moving puncture 
BH. We summarize our Bondi simulations in Table Hill 




t/Ji 

FIG. 13: Position of the BH puncture versus time for our 2-1-1 
axisymmetric runs with numerical resolutions of M/32 (run 
MP2a; solid), X/24 (run MP2b; dotted), M/16 (run MP2c; 
dashed), and M/8 (run MP2d; dot-dashed). In the bottom 
panel, we show the difference in position versus time for the 
same runs, noting that we see second order convergence in 
the position of the puncture. 



TABLE III: Summary of our matter evolution simulations. 
Key equations for the static Bondi solution, denoted by "SB" , 
are summarized in AppendixlB] Moving Bondi solutions, con- 
structed as described in Appendix lB 31 are "boosted" to have 
a linear velocity = 0.1 



Name 


Initial Data 


Grid 


Resolution 


tp 


2+1 w/Fisheye; \x}i^^\ < 15M, \x'\ < 136M 


SB2a 
SB2b 
SB2c 


Static Bondi 
Static Bondi 
Static Bondi 


360 X 360 
240 X 240 
120 X 120 


M/24: 

M/16 
M/8 


lOOX 
1007W 
1007W 


3-1-1 w/Fisheye; 


x'fishl < WM, \x'\< 56M 


SB3a 
SB3b 
SB3c 
SB3d 


Static Bondi 
Static Bondi 
Static Bondi 
Static Bondi 


240' X 120 
200' X 100 
160' X 80 
120' X 60 


M/12 
M/10 
M/8 
M/6 


lOOX 
lOOM 
lOOM 
lOOA^ 


3-fl w/Fisheye; = 0.1, \x}i,,, \ < lOM, 


< 56X 


MB3a 
MB3b 
MB3c 
MB3d 


Moving Bondi 
Moving Bondi 
Moving Bondi 
Moving Bondi 


240' X 120 
200' X 100 
160' X 80 
120' X 60 


M/12 
M/10 
M/8 
M/6 


lOOX 
lOOM 
lOOM 
lOOM 



A. Stationary Bondi solutions 

The analytic solution describing Bondi accretion onto 
a stationary BH is usually given in Schwarzschild coordi- 
nates (see, e.g., Appendix G of [79] for a detailed descrip- 
tion, and Appendix IB II for a summary) . To construct 



initial data for our dynamical simulations, we transform 
this solution into isotropic coordinates, as described in 
Appendix lB 21 Since isotropic coordinates become singu- 
lar on the horizon ai r = M /2, so does the fluid velocity 
when expressed in these coordinates. We therefore adjust 
the fluid initial data artificially in the immediate vicinity 
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of the BH. Specifically, for the rest-mass density, we fit a 
quadratic function between r — and r = so that 
its radial derivative matches the analytic Bondi solution 
at r = M. and the derivative goes to zero at r = M./2. 
We note that for stationary BHs, M — AIadm- Inside 
the horizon at r = A4/2, we set the density equal to a 
small positive value at the origin, plus a term with radial 
dependence oc 1— cos(27rr/A^) used to establish a smooth 
fit. For the velocity, we simply set u = u\^^j^ x [r/M) 
for r < M. Since the flow is supersonic and directed 
inward, this does not affect the exterior solution, and 
even within r < A4 the fluid solution quickly settles into 
an equilibrium flow as our spacetime slicing evolves to- 
wards the late-time solution that penetrates the horizon 
smoothly. 

To match to the stationary Bondi solution, in which 
the self-gravity of the gas is negligible, we require that 
the mass accretion rate multiplied by the integration time 
- which we choose to be if = lOOA^ - remain small with 
respect to the mass of the BH. Thus, we set M — 10^^ 
for all runs shown here. We set the sonic areal radius 
to Vs ~ 10 Ai. The proper infall time required for the 
fluid to travel from r « 9A4 (r^ in isotropic radii) to the 
horizon a,t rh = M/2 is 23 Ai, so that we evolve for just 
over four freefall times. The gas is at rest at infinity, with 
a uniform density of poo = 6.2 x 10~®. The polytropic 
index is chosen to be n = 3 (thus, the adiabatic index is 

r = 4/3). 

For the gauge conditions used in moving puncture sim- 
ulations, the Bondi solution is not time-independent. 
Similarly to the vacuum solutions described in Sec- 
tion |TTT1 the evolution passes through a transient, time- 
dependent phase, and then settles down into a new equi- 
librium. This new equilibrium solution describes the 
usual Bondi solution, but expressed in a different co- 
ordinate system. To analyze this solution and com- 
pare it with the analytical solution (which is given in 
Schwarzschild coordinates), we therefore need to com- 
pare invariants. 

One such invariant is the rate of change of the fluid 
rest density po a-s measured by an observer moving with 
the fluid, 

po = dpo/dr = uf^d^po, (19) 

where is the fluid's 4-velocity. 

In Fig. 1141 we show dpo/dr measured along the z-axis 
(perpendicular to the symmetry axis) as a function of 
the rest-mass density itself. In the top panel, we show 
the results for our highest resolution 3+1 simulation, run 
SB3a with A^/12 spacing, whereas in the bottom panel 
we plot values for run SB2a (A^/24). In both cases we 
plot the exact solution as points, along with our numer- 
ical profiles at t = 20M, t = 60M, and t = lOOM. We 
find good agreement throughout the evolution with vari- 
ations of no more than 10% and 4% respectively for the 
3-1-1 and 2-1-1 simulations. 

As an additional test we compute the average areal 
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FIG. 14: Time rate of change of the rest-mass density, dpa / dr, 
as a function of the rest-mass density. We show results for our 
highest resolution 3-1-1 (top panel) and 2+1 (bottom panel) 
calculations. The exact solution is represented by points, 
along with our profiles at t = 20M (dotted curve), t = 60M 
(dashed), and t = lOOX (sohd). 

radius of isodensity surfaces, defined by 

in terms of the surface's proper area A. Evidently, this is 
again a coordinate-independent quantity. For the essen- 
tially spherically symmetric isodensity surfaces consid- 
ered in this Section the average radius must of course be 
equal to its local value, but the definition Eq. I|20p gener- 
alizes to the non-spherical configurations in Section llVBI 
In Fig. [151 we show the average radii of the same iso- 
density surfaces described in Fig. [TH following the same 
conventions. We note that it is possible to spot the phase 
transition between subsonic and supersonic behavior at 
the sonic radius rs, marked by an arrow, which causes 
a shift in the power-law index of the radius-density re- 
lation. We also marked the horizon at Vh with a second 
arrow. Again, we see that the results are stable for a 
long period, with variations of no more than 5% and 3%, 
respectively. 

To test the convergence of our implementation of rel- 
ativistic hydrodynamics we have performed Bondi evo- 
lution calculations using the same set of numerical res- 
olutions used previously for the vacuum puncture calcu- 
lations described in Sec. IIIIBI In Fig. [121 we show ta 
for runs with differing numerical resolutions at the same 
fixed density value, chosen to be po — 3.8 x 10"'', which 
lies near the center of our logarithmic range and cor- 
responds to a location close to the sonic radius. We see 
that matter variables converge to second order in the grid 
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FIG. 15: Rest-mass density as a function of the areal radius of 
the corresponding isodensity surface, shown on a log-log scale. 
We show both our highest resolution 3-1-1 run (run SB3a; top 
panel) and 2+1 run (run SB2a; bottom panel). The density 
values at which we compute the surfaces are the same as in 
Fig. 1141 as are all other conventions, though we note that the 
density axis have been flipped. We mark the BH horizon 
and the sonic point of the flow rs with arrows, noting that we 
can see evidence for the well-known phase transition in the 
density-radius relation at the sonic point. 



resolution, just as the field variables do. In Fig. [TH we 
note that the extrapolated solution does seem to expand 
slowly over time, but that this effect represents approx- 
imately a 1% change over a period of t = 100A4. This 
effect is caused by the presence of the outer boundary. At 
sufficiently early times, our numerical solution converges 
to the analytical solution in regions that are causally dis- 
connected from the outer boundary. 

As we discussed in Section III C\ the point r — cor- 
responds to a surface of finite, positive areal radius (see 
(40l . l4l|). and represents a coordinate singularity only. 
Matter reaching this point therefore represents matter 
crossing a surface of a certain finite areal radius inside 
the horizon. Accordingly, all fluid quantities should re- 
main finite at r = 0. In the top panel of Fig. [T71 we 
show the rest-mass density as a function of areal radius 
for our highest-resolution 3+1 stationary Bondi result 
at different times, along with the analytical solution. 
We eliminate from the figure the three innermost grid 
points, where the hydrodynamical "fix" we apply (recall 
Sec. IIIEp directly affects the values of the primitive hy- 
drodynamical variables computed from the conserved set. 
We see that the flow extends smoothly within the hori- 
zon, extending inward to nearly the asymptotic limiting 
value of TA = 1.31A^ [lO]. Thus, the matter maintains a 
regular flow pattern into the BH, as we would expect, re- 
maining finite and well-behaved indefinitely. In the bot- 



FIG. 16: The average areal radius of an isodensity surface 
with density po ~ 3-8 x 10~^ as a function of time for 3-1-1 
runs with different numerical resolutions: runs SB3a {M/12; 
solid), SB3b {M/10; dotted), SB3c {M/8; dashed), and SB3d 
{A4/6; dot-dashed). In the bottom panel, we show the pair- 
wise differences between runs. 

tom panel of the figure, we show results at t ^ lOOA^ 
for runs of varying resolution, showing that we converge 
toward the analytic solution as we increase our numer- 
ical resolution, and that the physical region affected by 
the hydrodynamical fixes decreases in size as the resolu- 
tion increases. We find the convergence is second-order 
at larger radii, and approximately first-order nearer the 
puncture where differencing errors across the puncture 
and the hydrodynamical fixes impose small-scale oscilla- 
tions in the density. 



B. Bondi solutions for moving punctures 

We would also like to study our code's ability to han- 
dle matter in the presence of a moving puncture. For 
these purposes we study the Bondi solution as viewed 
in a frame in which the Schwarzschild BH puncture is 
moving. Our method is described in Appendix IB 31 

In this Section we consider BHs with a velocity of 
^ / J\/[ = 0.1, and let the BH start at a coordi- 
nate location of a; = — 2.5A^. As in the stationary Bondi 
case discussed in Sec. IIV Al above, we evolve our 3-f 1 
calculations for a duration t — lOOA^, equivalent to ap- 
proximately 4 sonic radius freefall times, at which point 
the BH has moved to a coordinate location oi x — 1.8A1. 
To visualize the resulting evolution, we show density con- 
tours with overlaid arrows representing the velocity field 
at t = 50A^ and 1007W for run MB3a in Fig. [HI We 
see a clear pattern of translation as the entire solution 
evolves. 
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FIG. 17: The top panel shows the rest-mass density po as a 
function of areal radius ta for run SB3a at times t = 50A4 
(dotted), t = 75M (dashed), and t = lOOX (dot-dashed), 
as well as the exact Bondi solution (solid). We exclude the 
three innermost gridpoints from each run, which are directly 
affected by our hydrodynamical "fixes". We see the solution 
remains smooth and accurate across the horizon, at ~ 
rh = 2M, which is marked with an arrow, remaining finite 
everywhere and nearly constant in time. In the bottom panel 
of the figure we show results at t = lOOA^ for runs with 
three different numerical resolutions, showing the expected 
convergence toward the exact solution. 

In Fig. [TO] we show as in Eq. (pO)) . and dpo/dr 
from Eq. (|19p as a function of density po for simulation 
MB3a, our highest resolution moving Bondi run (com- 
pare Figs. [HI and [T51 for stationary Bondi solutions). As 
before, the exact solutions are given by the points. We 
find that ta agrees with the analytical solution to within 
about 3%. The co- moving time derivative of the density 
dpo/dr, on the other hand, shows larger deviations of up 
to 20%. In Fig. [201 we show the convergence behavior 
in the average areal radius of the isodensity surfaces as a 
function of time for our runs, following the same conven- 
tions as Fig. [ini We see the same convergent behavior as 
in the vacuum and stationary Bondi cases: second-order 
convergence followed by the appearance of higher-order 
correction terms at approximately t = 40A1 . 




FIG. 18: Contours of the rest-mass density for our mov- 
ing Bondi evolution run MB3a, shown as slices through the 
equatorial plane at f = 50 (top panel) and t = lOOTVf (bot- 
tom panel). Our simulation extends to cover the negative 
y-plane as well, but is visually indistinguishable from being 
completely symmetric. Density contours begin at po = 10~^, 
and are spaced logarithmically, 16 per decade. A velocity 
vector representing a magnitude v — 0.2c is shown above the 
figure for reference. 
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V. DISCUSSION AND FUTURE 
CALCULATIONS 

In this paper we have performed several test simula- 
tions involving the modeling of BHs within the moving 
puncture approach, both in vacuum and in the presence 
of a relativistic fluid. 

Our vacuum tests focus on evolutions of both station- 



FIG. 19: Average proper radius of isodensity surfaces (top 
panel) and time rate of change of the density along the z-axis, 
dpo/dr, (bottom panel) shown as a function of the rest mass 
density on a log-log scale for run MB3a. The exact solutions 
are shown as square points. Results are shown a.t t = 20 
(dotted), t = &QM (dashed), and t = IQQM (solid). 



ary and moving BHs. We find that the code is second- 
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FIG. 20: The average areal radius va of the isodensity sur- 
face with density po — 3.8 x 10"*^ as a function of time 
(top panel) for runs with different numerical resolutions: runs 
MB3a {M/12; sohd), MB3b {M/10; dotted), MB3c (X/8; 
dashed), and MB3d {A4/6; dot-dashed). In the bottom panel, 
we show the scaled pairwise differences between the runs. 



efforts to model relativistic BHNS binaries. We have 
previously constructed BHNS initial data, using the 
conformal thin-sandwich formalism and imposing quasi- 
equilibrium boundary conditions on the BH's horizon 
(see [1^ H^)- We have also performed preliminary dy- 
namical simulations in conformal gravitation (see |47l . 
lisj). We now plan to adopt the puncture method to 
simulate BHNS binaries fully self-consistently, similar to 
the calculations of [iBI, [46[ . 

In contrast to [i^, [4^, we plan to evolve quasi- 
equilibrium conformal thin-sandwich initial data. This 
poses the conceptional problem that the initial data use 
excision to model the BH, while the dynamical evolu- 
tion requires initial data everywhere. As shown above, 
however, our experiments with single BHs demonstrate 
that we can replace the initial data in the BH interior 
with arbitrary sufficiently smooth functions without sig- 
nificantly affecting the late-time solution, or the evolu- 
tion, in the BH exterior. In fact, this suggests that mov- 
ing puncture simulations of BHBH binaries may also use 
initial data constructed in the conformal thin-sandwich 
formahsm (see, e.g., [13, [Ell), which are beheved to rep- 
resent true quasi-equilibrium configurations more accu- 
rately [5^. The evolution of conformal thin-sandwich 
binary BH initial data will be considered in detail in a 
forthcoming paper [60j . 



order convergent, as expected, but is limited primarily by 
the maximum numerical resolution we can achieve with 
our current second-order formulation (cf. the ADM mass 
convergence demonstrated in (soj at finer numerical reso- 
lution). To remedy this issue, we have introduced fourth- 
order spatial differencing into our code, results of which 
will be reported in [13] ■ We demonstrate that we can 
reproduce numerically the analytical solution of [i^l for 
an isolated stationary BH, in line with previous studies 
of stationary isolated punctures (see, e.g., (23l. [i^ 
We also demonstrate that we can artificially modify the 
initial data in the BH interior, and even violate the con- 
straints there, and still settle down to the same late-time 
asymptotic solution without significantly affecting the 
evolution in the BH exterior. While this result is not 
surprising, it is reassuring, and also has some implica- 
tions for future simulations, as we discuss below. 

To test the ability of the moving puncture approach 
to handle the flow of matter onto BHs, we study rela- 
tivistic Bondi flow both for stationary and moving BHs 
and again compare with the analytic solution. Our find- 
ings demonstrate that all fluid variables remain regular 
throughout and do not need excision. This result can 
be understood in terms of the studies of [H, [4l[, who 
demonstrate that the puncture represents a limiting sur- 
face of finite areal radius, and hence a coordinate singu- 
larity only. Moving puncture evolutions never encounter 
the BH's central spacetime singularity, and cover regular 
regions of the spacetime only. Hence all fluid invariants 
are regular throughout the puncture BH interior. 

This paper represents a stepping-stone in our group's 
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APPENDIX A: FISHEYE COORDINATES 

Fisheye coordinates are defined through a purely radial 
coordinate transformation, in which we define a fisheye 
radius r in terms of the physical radius r according to 

f - a r + 'S^ (a,-i - apSj / cosh((7- + Ri)/s.i) \ 
r-a^r ^ 2 tSinh{R, / s^) Uosh((r - i?,;)/s») J ' 

_ (Al) 
(see Eqs. (3) and (4) of [iMl)- Here the a,; coefficients 
determine the "stretching" of the radial coordinate in 
several different regions, labeled by the i's, the i?i's define 
the center of transition zones between these regions, and 
the Si's determine the width of these transition zones. 
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The derivative of this expression is given by 



We then have 



dr " tanh(i?i/si) 

/ tanh((r + Ri)/si) — tanh((r — Ri)/si) 



(A2) 



V 2 

which helps to understand how this transformation 
works. Assume that we arrange the Ri terms in such 
a way that neighboring transition zones do not overlap, 
i.e. Ri — Si > i?i_i+Si_i. Outside of the transition zones, 
i.e. for radii Rm-i + Sm-i < r < Rm ~ s.,„, the derivative 
df /dr is then given approximately by that region's coeffi- 
cient a,„_i. This is because the last term in parentheses 
takes a value approximately equal to zero for alH < m — 1 
and approximately one for i > m. The fisheye transitions 
act in many ways like an effective fixed-mesh refinement, 
with spherical transitions separating regions with differ- 
ent resolutions. Angles with respect to the origin remain 
unchanged by the transformation, and spheres centered 
on the origin transform into spheres, albeit with a differ- 
ent radius. We always apply the fisheye transformation 
in terms of the origin of our coordinates, regardless of 
the position of the BH in cases where it is moving across 
the grid. 

The fisheye transformation is purely radial. For Carte- 
sian coordinates we define 



X — X 



(;) 



(A3) 



The Jacobian of the coordinate transformation and its 
inverse are given by 



Ox^ f T\ 



dx 
dx^ 

dx^ 



x^x^ f dr r 

f'^ \df r 

x'^x^ f dr f 

dr r 



(A4) 
(A5) 



To construct initial data on a fisheye grid, the most 
direct method is to evaluate all quantities at the physical 
coordinates represented by the point, followed by an ap- 
propriate coordinate conversion. Since the transforma- 
tion is purely spatial, all time-components of 4-vectors 
remain unaffected, and we can restrict the transforma- 
tion to spatial components only, e.g. 



-i i - 



dx^ 
dx^ 



and 



dx' dx"^ 
dx^ dxi 



■ll-n 



(A6) 



(A7) 



The rest-density, for example, is the time-component of 
the fluid's density 4-vector, and is hence invariant under 
these transformations. 

The determinant of the Jacobian is given by 



J 



det{dx'/dx^) = (^) 



r\'^ dr 
df 



(A8) 



7 = det(7j) = J^7, 



(A9) 



indicating that the determinant of the metric is a tensor 
density of weight 2. Using a conformal transformation 
7ij = ip^Jij so that the determinant of the conformally 
related metric is unity, 7=1, implies 



2„/,12 



(AlO) 



meaning that the conformal factor is a tensor density 
of weight 1/6. The conformal factor and its logarithm 
(j) = Inip then transform according to 



$ = \n^ = iln f^W iln ( ^ ) -l-lnV-. (A12) 
3 Vr/ 6 \dr J 

Our conformal field quantities, "fij and Aij, both trans- 
form according to the relation. 



r\4/3 fdrV'^ dx' dx"" 



dr J dx^ dx^ 



■llr, 



(A13) 



The most complicated transformation is that of f% given 

by 



p ^ j2/3^fl _ ~l' 

dx' 



IdJ^dx^ 2/3 5^^^ 

2 dx'' 'dx''^ 



dx'dx"- 
,(Ai4) 

A time derivative of this term appears in the shift evo- 
lution equation (jllbp . Instead of evaluating this term 
exactly, we found it convenient to replace this equation 
with 



dtB' 



dtt' - vB\ 



(A15) 



since calculating the time derivative of the quantity 
r* requires taking a complicated spatially-varying lin- 
ear combination of time derivatives of the spatial met- 
ric. Within the different fisheye regions (i.e. outside the 
transition zones), our expression reproduces the "non- 
shifting-shift" condition, and it is equivalent to Eq. (jllbp 
in the area surrounding the BH itself. Elsewhere this 
modification affects only the coordinates, and not any 
physical quantities. 

Due to the asymptotic behavior of various quantities in 
fisheye coordinates, we must modify our boundary con- 
ditions in some cases in order to reproduce the desired 
behavior in physical coordinates. In all cases where out- 
going wavelike boundary conditions are used, we evaluate 
all radii in physical coordinates, not fisheye coordinates. 
The conformal factor which does not asymptotically 
approach unity in fisheye coordinates, is converted into 
the physical coordinate expression using Eq. (jA12p . at 
which point boundary conditions are applied and the ex- 
pression is converted back into fisheye. 
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APPENDIX B: 



THE RELATIVISTIC BONDI 
SOLUTION 



A thorough derivation of the exact analytic relativistic 
Bondi solution may be found in Appendix G of |79l] . Here 
we briefly review some of the basic features and most 
relevant equations. We assume that a BH of mass M is 
placed within an infinite cloud of gas that has a rest-mass 
density poo and fluid 3-velocity = at spatial inflnity, 
r — > oo. We take the gas to be adiabatic with adiabatic 
index V. We can then solve the equations of relativistic 
hydrodynamics to find the stationary spherical accretion 
flow onto the black hole. 



1. Review of key equations 

For convenience, we will derive the equations deter- 
mining the flow in Schwarzschild coordinates, and then 
convert these to the isotropic coordinates used through- 
out this paper. We denote Schwarzschild radii f, the 
4- velocity m", and the inwardly directed radial compo- 
nent of the 4- velocity u = —u''. In terms of these we can 
recast the relativistic continuity and Euler equations in 
conserved form 



1 - 



2M 



M — const. 



= h^, ~ const. 



(Bl) 
(B2) 



Here we define the specific enthalpy h = \ + e + P/ 
where e is the internal energy of the fluid. We will assume 
a gamma-law equation of state 

P = (r - l)poe, (B3) 

for which the enthalpy is h = 1 + Fe. Inserting the latter 
into (HH yields 

^=(r-l)po^. (B4) 

For adiabatic flow the gamma-law equation of state (|B3[) 
implies the polytropic relation P — kPq, where k is a 
constant. Combining this with (|B4p yields 



i^Po (r- l)po 



1 



(B5) 



For the gamma-law equation of state (EOS) (IB3|) the 
speed of sound is given by 



1 



dP 



1/2 



TP 

poh 



1/2 



(B6) 



Combining the above expressions we then find the fol- 
lowing relations between the enthalpy and the speed of 
sound 

a2 = (F-l)^, (B7) 



h 

,2 



1 



F - 1 - a2 



(B8) 



All smooth solutions to the conservation laws (|B1[) and 
(|B2p satisfying the EOS (|B5p must pass through a sonic 
point, since the flow is subsonic at large radii but must 
be supersonic at the horizon. It can be shown that at the 
sonic point the radial velocity its must satisfy 



.2 M 
2rs 

and that the speed of sound at the sonic point is 

,-•,2 



1 - 3u2 



(B9) 



(BIO) 



The accretion rate for the transonic solution is given 
uniquely by 

M = 47rp,u,f2 = 47rA,MVooai^ (BU) 

where poo and Ooo are the asymptotic density and sound 
speed, respectively, and As = As(F) is tabulated in Ta- 
ble 14.1 of [ll for values 1 < F < 5/3. 

In Fig. [211 we show the particular Bondi solution we 
use throughout this paper, corresponding to a F = 4/3 
EOS with a transonic flow satisfying M = 10^'' and 
fs = lOM (for convenience, we set M = 1). In terms 
of these parameters, the polytropic constant is given 
by K = 7.56, and the asymptotic rest-mass density by 
Poo — 6.2 X 10~^. In the top panel, we show the rest-mass 
density po as a function of both Schwarzschild radius f 
(dashed curve) and isotropic radius r (solid curve). Note 
that the isotropic solution terminates at = 0.5M, since 
the interior Schwarzschild solution is mapped through the 
throat of the BH onto the other sheet of the topology. In 
the second panel, we show the value of m° as a function of 
the two coordinate radii, showing the divergence at the 
horizon. In the third panel, we show the radial compo- 
nent of the respective 4- velocities, u{r) and u(f). Finally, 
in the bottom flgure, we show the radial component of 
the 3-velocities, v = \v^\ = {u^l/u^, seeing that in both 
cases this quantity goes to zero at the horizon, since the 
lapse vanishes there. We note, for clarity, that the rest- 
mass density po and u remain finite and smooth through 
the horizon in the Schwarzschild metric, becoming sin- 
gular only at the physical singularity at the origin. On 
the other hand, because of the coordinate singularities 
present in the Schwarzschild metric, the time-component 
of the 4- velocity u° and the radial 3-velocity both di- 
verge at the horizon. 



2. Transformation to isotropic coordinates 

To convert the Bondi solution from Schwarzschild to 
isotropic coordinates, we only need to transform the 
radii and velocities. The rest-mass density - as a time- 
component of the density 4-vector - is invariant under 
purely spatial coordinate transformations. The radial 
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FIG. 21: The Bondi solution for matter accreting onto a BH, 
where we set the fluid EOS to be a F — 4/3 polytrope, 
the sonic radius as fs = lOM, and the mass accretion rate 
M = 10~*. In the top two panels, we show the rest-mass den- 
sity and ii". Solid curves show the quantity as a function of 
isotropic radii, whereas dashed show Schwarzschild radii, with 
the transformation given by Eq. (|B13|I . In the third panel, we 
show the radial component of the 4-velocity in both coordi- 
nate systems, with the transformation law given by Eq. HB14|) . 
In the bottom panel, we show the radial component of the 3- 
velocities for both coordinate systems. We note that ^ oo 
and — > in either coordinate system as we approach the 
BH horizon, located at Vh ~ 0.5M in isotropic coordinates 
and fh = 2M is Schwarzschild. 



transformation is given by 



r = r[l + -) , 



f-M + ^r{r - 2 A/) 



Since the flow is purely radial, we have 



(B12) 
(B13) 



u = —u = —u — = u 

af \ dr 



(l-M/2r)(l + M/2r)' 



(B14) 



The 4-velocity's time-component u*^ remains unchanged 
under the transformation, which means that the radial 
component of the 3- velocity — /v? transforms in 
the same way as the that of the 4-vclocity. We find the 
3- velocity from the normalization 



-1. 



3. Moving Bondi solutions 



(B15) 



One approach to construct a moving Bondi solution is 
to view the solution of the previous section from a frame 
in which the BH is moving and described by a moving 
BH puncture spacetime. Our strategy is to adopt a mov- 
ing puncture spacetime and regular matter density and 
velocity profiles, with correct Bondi fiow outer bound- 
ary conditions. We then allow the solution to come to 
steady-state and use invariant flow variables to compare 
with the stationary Bondi flow solution. 

Specifically, we take the initial metric coefficients to 
be the vacuum moving puncture solution, as described 
in Sec. IIIBI We take the initial density at any coordi- 
nate point to be approximately the stationary isotropic 
Bondi solution. Finally, we compute the initial velocity 
field using the special-relativistic transformation law for 
velocities. Denoting the stationary Bondi solution with 
V, and assuming a "boost" velocity Vh 
direction, we have 



vf in the x 



V 

c 
c 
c 



- c "^vfV^) ' 



-2yxyxy 



(B16) 
(B17) 
(B18) 



where 7f, = (1 — w^/c^)^^/^ and the local value of the 
speed of light is c = a/ip"^. Here, Vb = P/M, where P is 
the momentum of the puncture. 

The resulting initial data matches the Bondi solution 
only approximately. However, all deviations propagate 
away quickly, leaving behind stationary Bondi flow onto 
a BH. 
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